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Results from an energy loss model that includes thermal fluctuations in the energy loss for heavy 
quarks in a strongly-coupled plasma are shown to be qualitatively consistent with single particle 
data from both RHIC and LHC. The model used is the first to properly include the fluctuations 
in heavy quark energy loss as derived in string theory and that do not obey the usual fluctuation- 
dissipation relations. These fluctuations are crucial for simultaneously describing both RHIC and 
LHC data; leading order drag results without fluctuations are falsified by current data. Including 
the fluctuations is non-trivial and relies on the Wong-Zakai theorem to fix the numerical Langevin 
implementation. The fluctuations lead to surprising results: B meson anisotropy is similar to that for 
D mesons at LHC, and the double ratio of D to B meson nuclear modification factors approaches 
unity more rapidly than even predictions from perturbative energy loss models. It is clear that 
future work in improving heavy quark energy loss calculations in AdS/CFT to include all energy 
loss channels is needed and warranted. 
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I. INTRODUCTION 

The goal of heavy ion collision research is to map out 
the many-body properties of quantum chromodynamics 
at extreme temperatures. Simultaneously, this research 
provides insight into the confinement transition of QCD 
matter and yields information on the earliest moments of 
the Universe’s existence: currently, no other experimen¬ 
tal program can reach back as far in time [1]. 

The fireballs created in central RHIC and LHC heavy 
ion collisions reach temperatures of a trillion degrees [2- 
5]; at this temperature the many-body dynamics of col¬ 
ored material is very different from that of normal nu¬ 
clear matter [6]. The best way of theoretically under¬ 
standing this material is not yet known. Simple ques¬ 
tions such as “What are the relevant degrees of freedom 
for the system?” do not yet have a clear answer. At the 
most basic level, it is not clear how to reconcile exper¬ 
imental results that, naively, simultaneously suggest a 
strongly-coupled plasma that evolves hydrodynamically 
and a weakly-coupled gas of slightly modified quarks and 
gluons that yields a plasma relatively transparent to hard 
probes. 

In order to make progress as a held it is necessary to 
attempt to effect such a reconciliation of ideas of the na¬ 
ture of the plasma into a consistent theoretical picture 
that yields quantitatively consistent comparisons with 
experimental data. Focusing on two major avenues of 
investigations, low transverse momentum (low-pr) ob¬ 
servables that are described well by near-ideal relativis¬ 
tic hydrodynamics switched on very soon after the ini¬ 
tial collision overlap [2-5] and high-p-p observables as¬ 
sociated with hard production described by perturba¬ 
tive quantum chromodynamics (pQCD) [7-9], the low-pT 


* wa.horowitz^Quct.ac.za 


observables appear readily understood within a strong¬ 
coupling paradigm for the plasma in which calculations 
are performed using the AdS/CFT correspondence [10- 
12] while the high-pT observables appear readily under¬ 
stood within a weak-coupling paradigm for the plasma in 
which calculations are performed using pQCD [7]. 

Specihcally, leading order derivations exploiting 
the AdS/CFT correspondence found that quark-gluon 
plasma (QGP)-like systems have thermodynamic prop¬ 
erties such as entropy density similar to those seen from 
the lattice [6, 10]; rapidly hydrodynamize after a heavy 
ion collision-like event [13], i.e. nearly ideal relativistic 
hydrodynamics approximates well the dynamics of the 
system soon after a heavy ion collision as is inferred from 
data; and have a very small shear viscosity-to-entropy ra¬ 
tio r]/s ~ 0.1 in natural units [10, 11], again as inferred 
from hydrodynamics studies of data [2-5]. At the same 
time, leading order derivations based on pQCD show dif¬ 
ferences to thermodynamic properties [14], do not hydro¬ 
dynamize rapidly [15], and have a viscosity-to-entropy 
ratio p/s ~ 1 [16] at least an order of magnitude larger 
than that inferred from hydrodynamics. 

On the other hand, energy loss models for light and 
heavy flavored particles based on calculations using lead¬ 
ing order pQCD show broad qualitative agreement with 
data from both RHIC and LHC [7-9]. At the same time, 
early light flavor [17] and leading order heavy flavor [18] 
energy loss models using AdS/CFT found or suggested 
a massive oversuppression of high-pr leading particles 
compared to data. 

Much work has been done to improve perturbative 
field theoretic treatments of soft, or low-pp, observables. 
For example there are high-order thermal field theory 
calculations of thermodynamic properties of the QGP 
[19] and recent work suggested that pQGD can yield 
rapid isotropization times of 0.2 — 1.0 fm [20]. The 
small viscosity-to-entropy ratio inferred from hydrody- 
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namics remains a challenge for perturbative calculations; 
no pQCD calculation has yet yielded ry/s ^ 0.1. 

An alternative research avenue is to improve the 
AdS/CFT treatment of high-pT observables. There are 
hybrid calculations in which, for instance, AdS/CFT is 
used to model the medium only [21-23] or only part of the 
energy loss [24], although in the latter case the coupling 
is taken to be quite small, A ^ 10“^ — 10“"^. Other 
work has attempted to apply AdS/CFT to the entire 
energy loss problem. A recent success [25] found that 
when the strong-coupling jet prescription is improved and 
the in-medium energy loss is renormalized, a reasonable 
value of A = 5.5 yielded a jet nuclear modification fac¬ 
tor R]^\{pt) quantitatively consistent with preliminary 
CMS data [26]. 

By allowing the variation of the intrinsic mass of the 
probe, heavy quark observables provide an extremely 
valuable additional test of energy loss calculations and, 
therefore, the properties of quark-gluon plasma. It has 
been very fruitful for the energy loss community to re¬ 
quire of energy loss models a simultaneous description of 
both light and heavy flavor observables. For example, 
[27] showed that pQCD-based energy loss models that 
only include the radiative energy loss channel could not 
simultaneously describe both suppression of particles as¬ 
sociated with both light and heavy flavors at RHIC. From 
the strong-coupling side of calculations, heavy quark ob¬ 
servables are especially useful as the theory of heavy 
quark energy loss is both more mature and under bet¬ 
ter control than that for light flavors [25, 28-32]. 

The work of this paper is to improve the energy loss 
modelling of heavy quarks under the assumption of a 
strongly-coupled medium strongly-coupled to the heavy 
quark. Previous calculations that used only the leading 
order strongly-coupled heavy quark energy loss, often re¬ 
ferred to as heavy quark drag, compared favorably to the 
measured nuclear modification factor Raa{pt) of elec¬ 
trons from heavy flavor decay at RHIC [33] but generally 
oversuppressed Raa{pt) for D mesons at LHC by a fac¬ 
tor of ~ 5 [18]. Another early comparison to RHIC data 
included fluctuations in the energy loss as dictated by 
the fluctuation-dissipation theorem [34]. These results 
were consistent with data from RHIC, but were never 
compared to data from LHC. A major deficiency of the 
approach of [34] is that the string theoretic derivation 
of the thermal momentum fluctuations for heavy quarks 
[ 21 ] demonstrated that the momentum diffusion is not re¬ 
lated to the momentum drag according to the fluctuation- 
dissipation theorem. Even worse, the fluctuations are 
non-Markovian; the fluctuations are due to colored noise, 
which means that there are non-trivial temporal correla¬ 
tions between the momentum kicks. 

In this paper the correct fluctuation spectrum is used 
and the non-Markovian nature of the momentum fluc¬ 
tuations is mitigated; it will be shown that by doing so 
one can qualitatively describe the data related to single 
heavy flavor observables at both RHIC and LHC simul¬ 
taneously. 


II. LANGEVIN ENERGY LOSS 


It was shown in [21] that the three-momentum p® of 
an on-shell heavy quark moving at constant velocity in a 
thermal bath evolves as: 

= -pPi + Ff" + F'[ , ( 1 ) 

where the drag coefficient [28-30] of a heavy quark of 
mass Mq in a plasma of temperature T with ’t Hooft 
coupling A is 


2Mq 


( 2 ) 


and the fluctuating momentum kicks are correlated as 


(E/'(ti)F/'(ti)) = KLP^PJg{t2 - h) (3) 

{F^iti)F]^(ti)) = KriS.j - p^pj)g{t2 - h), (4) 

where Pi = pi/\p\, 

kt = 7r-\/AT^7^^^ (5) 

KL = 7 ^/ ct , ( 6 ) 


and 5 is a function known only numerically. 

Since the fluctuations increase in magnitude with 
heavy quark velocity—in fact increase very, very quickly 
in the longitudinal direction (like 7 ^'^^), which is the 
most important direction for calculations of suppres¬ 
sion observables—one may naturally ask: at what speed 
should one expect the fluctuations to play a nontrivial 
role in heavy quark energy loss? This speed limit may 
be estimated by requiring that the momentum picked up 
via fluctuations over the time scale set by the drag coef¬ 
ficient be small in comparison to the total momentum of 
the heavy quark. One finds then an upper limit on the 
speed of the heavy quark given by 


Ml 

ry < = _ Q 

I I crtt 4X'^ * 


( 7 ) 


There is a well-known speed limit for the heavy quark 
drag setup in which the derivation of Eq. (1) was per¬ 
formed. A finite mass heavy quark is represented by a 
string connected to a D7 brane that partially fills the fifth 
dimension [35]. The mass of the heavy quark is related 
to the depth that the D7 brane fills the hfth dimension; 
the greater the depth, the lighter the quark. In a space 
with a black hole—dual to a plasma at hnite temperature 
[36] —the local speed of light decreases with the depth of 
the fifth dimension. Thus in this setup the endpoint of 
a string dual to a finite mass heavy quark that is con¬ 
nected to the bottom of the D7 brane has as a natural 
speed limit the local speed of light [ 21 ]; a heavy quark in 
the field may move faster than this speed, but then the 
setup is necessarily incorrect for this physical situation. 
Equivalently, the speed limit can be derived from self- 
consistency: in the setup in which the drag formulae are 
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derived, the heavy quark moves at a constant velocity. 
In order for the constant velocity to be maintained, one 
must provide power to the heavy quark to balance the 
momentum lost. If that power is provided by a gauge 
field, there is a maximum field strength before which 
Schwinger pair production begins; that maximum field 
strength limits the velocity with which the quark can be 
pulled through the plasma at a constant velocity [22]. 
Both lines of reasoning lead to exactly the same speed 
limit. This speed limit on the heavy quark setup leads 
to the requirement that 


7 < IcUt = (1 + 




4M2 
Ar2 ■ 


( 8 ) 


One may easily see that parametrically the critical 
speed is restricted most by the speed limit due to the 
requirement of constant heavy quark velocity, not due to 
the momentum kicks picked up from the medium; the 
former goes as while the latter goes as A°. How¬ 
ever, in the case of finite A ~ 0(10) one finds that the 
critical speed due to momentum fluctuations is actually 
smaller than that of the speed limit. It is therefore impor¬ 
tant to include these fluctuations in a phenomenological 
heavy quark energy loss model based on AdS/CFT that 
is compared to data, the goal of this work. 

The implementation of the fluctuations in Eq. (1) are 
nontrivial for two reasons: first, the noise is colored and 
second, the fluctuations are multiplicative. With respect 
to the former complication, if g was a Dirac delta func¬ 
tion, the momentum fluctuations in Eq. (1) would cor¬ 
respond to white noise; i.e. one momentum kick is not 
correlated in time with any other. However, it was 
shown in [21] that the noise is actually colored, in which 
case g is some nontrivial function of characteristic width 
tcorr ~ ■\JllT'■ With respect to the latter complication, 
multiplicative noise has correlations that depend on the 
dynamical variable, in this case pi. Unlike usual Rie- 
mann calculus, stochastic Ito integration is sensitive to 
the precise location at which the integrand is evaluated 
within a time step; different choices of location within 
the time step yield different results for the integral [37- 
40]. Note that these different results for the integral are 
truly fundamental and not due to finite size time steps. 
Unfortunately, the work of [21] did not specify at which 
place within the time step the derived fluctuations should 
be applied. Eig. 1 shows the significant differences in the 
momentum distributions of charm quarks in a thermal 
plasma of T = 400 MeV which is moving at n = 0.9 c 
in the 2 direction for three different but common choices 
made for the location within the time step to apply the 
momentum fluctuation. In particular the results for the 
ltd (pre-point, for which the integrand and the fluctua¬ 
tion are evaluated at the earliest time within a time step), 
Stratonovich (mid-point, for which the integrand and the 
fluctuation are evaluated at the middle of a time step), 
and Hanggi-Klimontovich (post-point, for which the inte¬ 
grand and the fluctuation are evaluated at the latest time 
within a time step) rules are shown. The differences be¬ 


tween rules are especially pronounced in the tails of the 
distribution (c.f. Fig. 1 (b)), which tends to be the most 
important for the purposes of energy loss modeling. Note 
that none of the distributions produced by the different 
integration prescriptions for the momentum fluctuations 
given by AdS/CFT are exactly thermal. The deviation 
away from the thermal distribution can be attributed at 
least in part to the modification of the usual dispersion 
relation of the heavy quark due to its coupling to the 
strongly-coupled thermal medium [29]. In Fig. 1 the rel¬ 
ativistic Maxwell-Boltzmann, or Jiittner, distribution is 
shown and compared to a distribution produced by a 
Langevin evaluation whose diffusion coefficient is related 
to its drag via the relativistic fluctuation-dissipation re¬ 
lation [41]. This comparison between the analytic and 
numerical results is a non-trivial cross check of both the 
analytics and numerics: the distribution along the direc¬ 
tion of plasma fluid flow. Fig. 1 (b), is not a trivial boost 
of the distribution along a direction for which the plasma 
is at rest. Fig. 1 (b); see Appendix A. 

Fortunately, both the colored and multiplicative com¬ 
plications of the noise can be handled. As for the first 
complication, in every physical situation, all noise is ac¬ 
tually colored (for example the air molecule that just 
collided with a floating pollen grain cannot immediately 
collide with the pollen grain again); the issue is whether 
or not the time scale over which the momentum kicks are 
correlated is small compared to the other relevant time 
scale(s) in the problem. For heavy quark energy loss as 
computed at strong coupling in AdS/CFT, the only other 
relevant time scale is determined by the drag coefficient 
p, from Eq. (I). One may check that 

1 T 

tcorrl^ ^ —'\/Ay^ < 1 (9) 

2 Mq 

so long as 7 < from Eq. (8). Thus we are able 

to approximate the colored noise as white noise so long 
as the heavy quark is moving slower than the speed limit 
imposed by the setup of the derivation. As for the second 
complication, there is a theorem due to Wong and Zakai 
[42] that proves that in the limit that the autocorrelation 
time for noise goes to zero, the correct stochastic integral 
is the Stratonovich one, which is to say that the integrand 
and noise terms are evaluated at the midpoint of the time 
step. As just noted, for the specific calculation of interest 
in this work, the autocorrelation time is small so long as 
the heavy quark moves more slowly than the speed limit 
'ycriti Eq. (8). 

In the energy loss model used in this paper, the 
stochastic differential equation Eq. (1) was solved nu¬ 
merically using the Euler-Maruyama scheme [40]. De¬ 
spite its simplicity, the Euler-Maruyama scheme con¬ 
verges strongly, and the scheme was sufficient for the 
purposes in this paper. (Higher order strongly convergent 
schemes exist, e.g. Milstein, but are highly non-trivial to 
implement for motion in more than one dimension.) In 
any implementation of a stochastic integral, whether Ito, 
Hanggi-Klimontovich, or anything in-between, one may 
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FIG. 1. (a) The p“^-differential distribution of rric = 1.5 

GeV/charm quarks strongly-coupled to a strongly-coupled 
plasma in three spatial dimensions at T = 400 MeV which 
is moving at 0.9 c in the a direction. The effect of the 
thermal drag and fluctuations on the heavy quarks from 
strong coupling is found from implementing Eq. (1) using 
the Ito (pre-point), Stratonovich (mid-point), and Hanggi- 
Klimontovich (post-point) stochastic integrals. An implemen¬ 
tation for which the fluctuation-dissipation theorem holds is 
also shown and compared to the result expected, the Jiittner 
distribution, (b) The same as (a) but for the p^-differential 
distribution. 


actually evaluate the integral in any other choice with 
the addition (or subtraction) of the appropriate terms 
[38, 40]. In particular an ltd (pre-point) stochastic dif¬ 
ferential equation 

dp^ = a*(t, p^{t))dt + (t, p^{t))dWj (10) 


where 


«'(*./(<)) =a‘(i./W) 

+ ( 12 ) 

For ease of numerics, all results shown in this work had 
the integrand and noise terms evaluated at the earliest 
moment in the time step (ltd, or, pre-point). Thus in 
the local fluid rest frame the momentum components of a 
heavy quark at the next time step, were found from 

the momenta at the current time step, p'^, by interpreting 
Eq. (1) as a Stratonovich stochastic differential equation, 
modified by Eq. (12) to be implemented as an ltd SDE 
in the Euler-Maruyama scheme with the result that 


/ 1 , 

Pn+l =(l-pdt'+ -Kdt'{ 

where p is given by Eq. (2); dt' is the time step dt boosted 
into the local rest frame of the fluid; dt' = dt/j] k = 
where T is the temperature of the fluid in its 
local rest frame; d is the number of spatial dimensions 
in the calculation (in this case we propagate the heavy 
quarks through backgrounds generated by VISHNU [3, 
4], which is a 2 -I- lU hydrodynamics code); the dWj are 
the uncorrelated, Gaussian Wiener kicks with mean 0 
and standard deviation 1; and 


= Vdt' Kyi/-* 


p^pO 


(7^ + 1) m ? 



(14) 


In a calculation that follows the trajectory of the heavy 
quark from production onwards, at each time step the 
code boosts into the local rest frame of the fluid, evalu¬ 
ates the change in momentum, boosts back to the lab 
frame, and propagates the heavy quark in coordinate 
space according to 


<+i = 


E, 


n+l 


(15) 


dt was taken to be 1/150 x pmax, where Pmax is the drag 
coefficient at the center of the fireball at the thermaliza- 
tion time, the largest drag coefficient for any individual 
collision. This value of dt was found through trial and 
error to consistently yield stable static thermal distribu¬ 
tions such as those shown in Fig. 1. Variations of dt by a 
factor of 2 in both directions yielded unchanged results 
in trials. The code was written in C-|—1-, compiled using 
Clang, and run on a Mac laptop. 


with Wiener process dW^ is equivalent to the HI- ENERGY LOSS MODEL 

Stratonovich stochastic differential equation of the form 

The results presented here are from a fully Monte- 
dp' = ad{t, p^{t))dt -\- {t, p^{t)) o dWj, (11) Carlo energy loss model; only the fragmentation from 
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heavy quarks into their decay products was handled 
through (numerical) integration of analytic expressions. 
In particular, the production spectrum for heavy quarks 
(from FONLL [43-45]), the spatial distribution of their 
production points (weighed by the Glauber binary dis¬ 
tribution [46]), and their initial direction of propagation 
(assumed uniform) were randomly sampled. Propaga¬ 
tion was through backgrounds generated by the VISHNU 
2-l-lD hydrodynamics code [3, 4], with momenta updated 
according to Eq. (1). (Pseudo-) random number genera¬ 
tion was performed using the Ran routine from Numerical 
Recipes [47]. Seed generation came from the routine de¬ 
scribed in [48] , which is appropriate for running multiple 
processes at or nearly at identical times on multiple pro¬ 
cessors^ . 

The production spectrum and fragmentation functions 
of the heavy quarks from FONLL were used [43-45] . The 
spectra were generated by the FONLL web page [49] for 
the relevant center of mass and rapidity ranges for the 
various experiments. The fragmentation functions from 
heavy quarks to heavy mesons were taken from [45]. For 
measurements of unspecified D mesons (and for electrons 
from D mesons), the FONLL approximation of D mesons 
being comprised of 70% and 30% 77+ mesons was used 
[49]. The semi-leptonic decay of heavy mesons to elec¬ 
trons is not given explicitly in any FONLL paper. The 
probability P{E) of a heavy meson decaying to an elec¬ 
tron of energy E in the rest frame of the heavy meson 
used in [50] was provided [51], but detailed instructions 
for implementation were not. One may calculate the dis¬ 
tribution of electrons in the lab frame from the semi- 
leptonic decay of a heavy meson; the details are given in 
Appendix B. In particular, Eqs. (B2) and (B7) - (BIO) 
are used throughout the rest of this work for the frag¬ 
mentation of heavy mesons to electrons and B mesons to 
J/ip mesons. 

The production spectrum of electrons originating from 
open heavy flavor \x\ p + p collisions at ^/s = 200 GeV 
from FONLL is compared to that measured by PHENIX 
[52] at RHIG in Fig. 2. One can see that at ultrarelativis- 
tic momenta, the FONLL predictions are consistent with 
the data within the experimental and theoretical uncer¬ 
tainties, which amount to approximately a factor of 2. 
At lower momenta the FONLL predictions deviate more, 
rising to a discrepancy of an order of magnitude for the 
smallest measured momentum bin of pr = 0.25 GeV/c. 
For electrons it appears that the production and/or the 
fragmentation description is not reliable below ^ 3 
GeV/c. (It seems likely that the issue is with the frag¬ 
mentation to electrons, with the discrepancy made much 
worse by the experimentally restricted rapidity range.) 

When computing the energy loss and fluctuations, one 


^ Random number generators are often seeded by the run time 
of the code; when multiple instances of the code are created at 
the same time, runs can have additional, unintended correlation 
from using the same seed. 



FIG. 2. (Top) A comparison of FONLL predictions [43- 
45] for the electron spectrum from open heavy flavor produc¬ 
tion in p -b p collisions at ^/s = 200 GeV to that measured 
by PHENIX [52]. (Bottom) Ratio of the PHENIX data to 
FONLL predictions; the solid black lines represent the un¬ 
certainty in the FONLL predictions due to scale variations. 
(Inset) Zoom in on the low-pr part of the ratio of data to 
FONLL predictions showing the full range of disagreement. 

needs a mapping from the QGD parameters to the M = 
4 SYM parameters. Results in this paper are shown using 
two prescriptions, which are denoted in the plots “as = 
0.3” and “A = 5.5,” in order to explore a reasonable 
breadth of possibilities (one hopes that a calculation in a 
theory dual to QGD would lie somewhere between these 
two prescriptions). For the former, the’t Hooft coupling 
is taken to be A = dTrOsW = drr x 0.3 x 3 and Tqcd = 
Tsym- For the latter, the prescription of [21] is followed, 
in which A = 5.5 and Tsym = TQcol'i^^'^■ The masses 
of the quarks were taken to be me = 1.5 GeV/c^ and 
mb = 4.75 GeV/c^. 

IV. RESULTS 

Fig. 3 shows the main results of this paper: predic¬ 
tions for RaaIpt) for electrons from open heavy flavor 
at RHIG, D meson suppression at LHG, and non-prompt 
J/'i/ suppression at LHG all compared to data [52-55]. 
When including only the leading order drag without fluc¬ 
tuations, the model significantly oversuppresses the nu¬ 
clear modification factor compared to data. But when 
the correct fluctuations are included, the predictions— 
















6 





FIG. 3. (a) Open heavy flavor electron Raa{pt) at RHIC 
measured by PHENIX [52] and STAR [53] compared to pre¬ 
dictions from heavy quark drag only (“LO”) and heavy quark 
drag with fluctuations (“NLO”) energy loss models, (b) D 
meson Raa{pt) at LHC measured by ALICE [54] compared 
to the strongly-coupled energy loss predictions, (c) Non¬ 
prompt J/ip meson Raa{pt) at LHC measured by CMS [55] 
compared to the energy loss model predictions. For curves 
labeled Qs = 0.3, A = drr x 0.3 x 3 and Tsym = Tqcd; for 
curves labeled A = 5.5, Tsym = ■ 


within their regime of applicability—are in qualitative 
agreement with all current suppression data. 

To expand on the point of the regime of applicability, 
while the electron prediction below pt ^ 3 GeV/c is 
not in particularly good agreement with data, the elec¬ 
tron production and/or fragmentation functions are un¬ 
trustworthy below this momentum scale. The regime of 
applicability is also crucial to note as it explains the dra¬ 
matic increase in Raa for electrons and D mesons for 
momenta beyond a few GeV/c (> 4 GeV/c for electrons 
and ^ 10 GeV/c for D mesons): self-consistently, one 
sees that at momenta of order the speed limit, the mo¬ 
mentum fluctuations begin to dominate at the momenta 
at which the setup of the problem is breaking down. The 
lack of inclusion of Larmor-like energy loss [56, 57] due 
to what would be significant deceleration becomes a poor 
approximation, as is the white noise assumption for the 
fluctuations^. When that important physics is included, 
the rapid rise in Raa{pt) above pt ^ 4 — 5 GeV/c will 
be tamed, although the magnitude of the effect is not yet 
known or estimable. 

In previous work [33], drag-only strong-coupling en¬ 
ergy loss gave a good quantitative description of RHIG 
electron data, but, when using the same parameters, pre¬ 
dicted an oversuppression of D mesons at LHG [18], due 
to the very sensitive p ^ nature of the drag. In this 
calculation there is already an oversuppression at RHIC. 
It is precisely this considerable sensitivity to temperature 
that leads to an oversuppression at RHIC in this calcu¬ 
lation: the VISHNU hydrodynamics [3, 4] is significantly 
hotter than the Glauber model used in the previous work 
[33]. 

One may also compare to other observables, for in¬ 
stance the azimuthal anisotropy of the D mesons, as 
shown in Fig. 4. The energy loss model including fluc¬ 
tuations qualitatively agrees with the experimental re¬ 
sults from ALICE [58], although with a generic under¬ 
prediction for the V 2 at intermediate-p-r values pr ^ 10 
GeV/c. This underprediction in D mesons is reminiscent 
of the underprediction of V 2 for light flavor observables 
in pQCD-based energy loss models [59, 60] and, like in 
the light flavor case, probably indicates a modification of 
the fragmentation process that enhances the anisotropy 
in this momentum range. Future measurements with in¬ 
creased statistics will, as in the light flavor case [17, 61], 
presumably show that the D meson V 2 decreases as a 
function of pt to levels similar to the energy loss predic¬ 
tions. Predictions for B meson V 2 from the energy loss 
model are also shown in Fig. 4; surprisingly, the model 
predicts B meson V 2 of a similar size as the D meson V 2 - 


^ Since the fluctuations increase dramatically with 7 , a rare few of 
the quarks can actually unphysically runaway to arbitrarily large 
momenta; the quarks fluctuate between increasingly large posi¬ 
tive and negative momenta for each time step. In the code, there 
is a cutoff of 1000 GeV/c, after which the quark is considered un¬ 
physical. The number of these “lost” quarks is 0(10“'*%). 
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FIG. 4. V 2 {pt) for D mesons as measured by ALICE [58] 
at LHC compared to the predictions from the heavy quark 
energy loss model including fluctuations. For comparison the 
fluctuating energy loss model predictions for B meson V 2 are 
also included. Curves labeled as = 0.3 use A = dvr x 0.3 x 3 
and Tsym = Tqcd', curves labeled A = 5.5 have Tsym = 
TqcdI^I^. 


It will be very interesting to see if measurements of the 
B meson vi are of a similar magnitude to the D meson 
V 2 as predicted here. 

Ideally, one would very much like some kind of smoking 
gun experimental measurement to distinguish between 
weakly-coupled energy loss mechanisms in the plasma 
compared to strongly-coupled energy loss mechanisms. 
The double ratio of the nuclear modification factor of D 
mesons to B mesons was proposed as just such a dis¬ 
tinguishing measurement [62]. It was shown in [62] that 
the leading order drag prediction for the double ratio is 
flat as a function of px and at a value of approximately 
niD/mB — 0.2. On the other hand, pQCD-based en¬ 
ergy loss models predict an insensitivity to the mass of 
the heavy quark when px ^ Mq, and thus the dou¬ 
ble ratio starts around 0.2 for px ^ 10 GeV/c and ap¬ 
proaches 1 from below by px ~ 100 GeV/c. It is nat¬ 
ural to ask whether the strong-coupling prediction for 
the double ratio oi D to B meson nuclear modification 
factors is affected by the inclusion of fluctuations. It is 
shown in Fig. 5 that fluctuations completely change the 
prediction for strong-coupling: the extremely rapid rise 
in Raa{pt) for D mesons due to fluctuations is reflected 
in an extremely rapid rise in the double ratio, with the 
double ratio actually going above 1 at p-r ~ 7 GeV/c. 
Of course it is precisely at this momentum range that the 
drag calculation is breaking down; an improved calcula¬ 
tion that includes the additional Larmor-like energy loss 
and autocorrelations in the fluctuations will bring down 
the double ratio prediction, although it is not yet clear 
by how much. It seems likely that, for single particle 
observables, the best experimental measurement to dis¬ 
tinguish between heavy quark drag-like energy loss and 
pQCD-like energy loss will come from B meson Raa{pt) 
alone: Fig. 3 (c) shows the drag (with fluctuations) pre¬ 



FIG. 5. The ratio of D meson Raa{pt) to B meson Raa{pt) 
from the strongly-coupled heavy quark energy loss model 
at leading order (dashed) and including fluctuations (solid). 
Curves labeled Qs = 0.3 use A = dvr x 0.3 x 3 and Tsym = 
Tqcd', curves labeled A = 5.5 have Tsym = TqcdIS^^*- 


diction decreases significantly with px and remains rela¬ 
tively constant while pQCD predicts that the B meson 
Raa{pt) will rise as a function of momenta for px ^ 25 
GeV/c. One might notice that there is no prediction for 
^aa''^'^ impart) given in this work. Since the strong¬ 
coupling energy loss model qualitatively describes the 
most central nuclear modification factor for non-prompt 
J/i/: mesons, and the model necessarily will predict that 
the nuclear modihcation factor goes to unity as collisions 
become more peripheral, it is difficult to imagine that 
the energy loss model will not follow the experimentally 
measured trend [63]. 


V. CONCLUSIONS AND OUTLOOK 

Hard probes, and, in particular, heavy quarks, pro¬ 
vide a crucial test of our theoretical understanding of 
the physics of quark-gluon plasma. Previous work [18] 
that included only the leading order energy loss [28- 
30] for heavy quarks strongly-coupled to a strongly- 
coupled plasma as calculated from the AdS/GFT corre¬ 
spondence showed significant, falsified disagreement with 
data [54, 64]. The results presented in this paper in Figs. 
3 and 4, for which energy loss fluctuations [21] are cor¬ 
rectly included for the first time in a realistic strong¬ 
coupling energy loss model for heavy quarks, are quali¬ 
tatively consistent with data from RHIC [52, 53] to LHC 
[54, 55, 58] within the regime of applicability of the en¬ 
ergy loss formulae, Eqs. (1) - (6). 

Along with the recent improvement in the treatment 
of light flavor energy loss in AdS/CFT [25], in which 
very good agreement between the suppression of jets 
from strong-coupling and the preliminary CMS data [26] 
was shown, strong-coupling can now explain within a 

































self-consistent picture of quark-gluon plasma dynamics 
a wide and impressive array of observables from low- 
momentum to high-momentum that includes the rapid 
applicability of hydrodynamics [13], the very small shear 
viscosity-to-entropy ratio [10, 11], the light flavor jet sup¬ 
pression [25], and the suppression of decay products of 
heavy quarks (this work). 

From the specific perspective of heavy flavor energy 
loss from AdS/CFT, the next steps include further test¬ 
ing theoretical predictions against data—especially in 
comparison with predictions from pQCD—and in im¬ 
proving the theory itself. 

Surprisingly, as shown in Fig. 5, the inclusion of fluc¬ 
tuations completely changes the leading order [62] qual¬ 
itative prediction of the double ratio of H to i? nuclear 
modification factors from the strong-coupling energy loss 
model: with fluctuations, the double ratio rises rapidly 
as a function of meson momentum, broaching unity at 
a momentum scale somewhat above where the calcula¬ 
tion is inapplicable due to the speed limit of the setup, 
Eq. (8). Although the rise is now more rapid than that 
predicted by pQCD [62], it is not known how much the 
momentum rise in the AdS/CFT prediction will be soft¬ 
ened as necessary corrections are included to extend the 
regime of applicability of the calculation to pt above a 
handful of GeV/c; it is therefore difficult at this stage to 
use an experimental measurement of the double ratio to 
test strong- versus weak-coupling physics. 

The anisotropy of heavy flavor decay products pro¬ 
vides a strong- versus weak-coupling energy loss compar¬ 
ison, although of limited utility. From strong-coupling, 
as shown in Fig. 4, B meson V 2 {pt) is predicted to be 
of a similar magnitude to D meson V 2 {pt)- Due to mass 
effects one expects pQCD to predict a smaller anisotropy 
for B mesons than for D mesons, but it does not seem 
that predictions for v^ipr) are extant in the literature. 
This prediction probably has limited use for distinguish¬ 
ing between strong- and weak-coupling dynamics as a 
precision differentiation of V 2 for D versus B mesons 
poses a significant experimental challenge. Addition¬ 
ally, it is not clear the extent to which non-perturbative 
hadronization effects differ for D versus B meson pro¬ 
duction in heavy ion collisions at relatively low < 10 
GeV/c momentum at which the AdS/CFT calculations 
are within their regime of applicability for D mesons. 

A likely much more useful prediction comes from the 
nuclear modification factor of B mesons. As can be seen 
in Fig. 3 (c), the strong-coupling energy loss model pre¬ 
dicts a large B meson suppression, Raa(pt) 0.1, out 
to pt 30 GeV/c; in contrast, pQCD predicts a much 
smaller B meson suppression with Raa{pt) 0.5 [18]. 
It is worth emphasizing that the strong-coupling energy 
loss derivations and implementations are under increas¬ 
ing theoretical control as the mass of the quark increases, 
so future measurements related to b quarks will be espe¬ 
cially interesting to see. 

Since momentum fluctuations play such an important 
role in strongly-coupled heavy quark energy loss, and 


these fluctuations grow characteristically with momen¬ 
tum, a comparison between theory and data related to 
correlations of pairs of heavy quarks would likely be very 
enlightening. The comparison would be most meaningful 
for the decay products of b quarks at high momentum, 
10 ^ Pt ^ 100 GeV/c: above ~ 10 GeV/c any effect 
that leads to the enhanced V 2 at intermediate-pT should 
be small, and one should be able to trust the theoretical 
calculation below ~ 100 GeV/c. The determination of 
these correlations from this energy loss model is left for 
future work. 

Improvements to the theoretical calculations that need 
to be done would increase the regime of applicability of 
the energy loss formulae and model, thus allowing for a 
meaningful comparison to more data. The heavy flavor 
electron predictions could be extended to lower momenta 
if their production in p -|- p collisions is brought under 
better phenomenological control. Both the heavy flavor 
electron and D meson predictions could be extended to 
higher momenta if the theoretical and numerical com¬ 
plications related to the speed limit, Eq. (8), could be 
overcome. It is likely that an entirely different approach 
to the energy loss calculation, both analytic and numer¬ 
ical, will be required: since the setup of a heavy quark 
moving at a constant speed is breaking down, presumably 
one must allow the heavy quark to slow down, which, in 
turn, will likely mean that full numerical solutions to the 
worldsheet evolution including thermal fluctuations will 
be needed. This important, non-trivial research is also 
left for future work. 
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Appendix A: The Jiittner Distribution 

The number of particles N of a relativistic gas in the 
rest frame of the gas in an n-spatial-dimensional box is 
given by the integral over the invariant phase space of 
the Jiittner distribution [41, 65, 66], the relativistic gen¬ 
eralization of the Maxwell-Boltzmann distribution, 

N = J d^xd^pAe-P°/^, (Al) 

in units of?i = c = fc=l, where p^ = \/ (p^ -|- ) is 

the relativistic energy of a particle of mass m, which can 
be 0. T is the well-defined temperature of the gas in its 
rest frame, and A is a normalization constant. 

Since we are interested in comparing to a numerical 
simulation with fixed number of particles N we may read- 
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ily determine A: 


where 


distribution of the JV particles is given by 

e-vVT 


f p”-idpe-Vp"+’”'/^, 

(A2) 

dN 

lo 


dPp 2 \2 tt'1') 


m^K n+i (m/T) 


(AS) 
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27r'^ 

iwy 


(A3) 


is the volume of a d-dimensional sphere; e.g., ^2 = 47r. 
The integral in Eq. (A2) may be computed analytically 
in terms of Gamma and modified Bessel functions, such 
that 

A=- —( —. (A4) 

2V\2ttTJ m-^K^im/T) ^ ^ 

Thus in the rest frame of the plasma the momentum 


Picking out one special direction, say the z-direction, one 
may find the p^-differential distribution of particles in the 
rest frame of the plasma: 


dN 

dp^ 


-n(—) 

2 \2ttTJ 
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[m/T) 

X J d"-V t e"V'PT+(p")"+™Vr_ 


If one is interested in the distribution of particles in a 
frame that is boosted along, say, the z-direction at a 
speed /? compared to the rest frame of the plasma, then 
one finds that in the primed coordinate system 


dN 

dp^' 
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(^Kr. (7y^(p^')2 -bTO^/r) 




, _^_ =Kr^ 

V (p^ )2 -b to2 2 


(A9) 


Appendix B: Semi-leptonic and J/ijj Fragmentation Functions 

The spectrum of electrons dN^/dp^{p^) (or J/ip mesons) produced in the lab frame by a spectrum of heavy mesons 
(in this case a D meson for specificity) that decays into an electron of energy E with probability P{E) in the rest 
frame of the heavy flavor meson is given by: 

= J dpT'^dE‘^:^{p!^)P{E)6(p^ - E^J-fl,{ sm{e) cos{(l)) + Pd) -b sin^(0) sin2(^!))) (Bl) 

poo pSma^ JA dN ^ 

= J^ dp? J sm{9)d6 J -^J{9, p)-^j^{p?)P{pD J{9, (t)))9{\sec{9)\J~^{9, p) - csch{yjnax)), 

where 


J{9, P) = 


jD(sin(9) cos(p) + Pd)^ -b sin^(d) sin^((/)) 



(B3) 
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Pmzu = max 


(O, ruD 


{Pt? - 


2 EjYiax 

^min — mO-X^O, (£^^0,3;)^ , (£773^73)^ 

^max — min(^7r/2, Re(sin ^{Pd + Pt/id Ermn))'j 
4^min — R'^(^crit (-^mm)) 

4*max — R'6(^/^crzt('^maa:)) (®^) 

= sin-^TPoE Pt/'IdE) 

(t)crit{E) = sin"^ ( -—^-— ( - 1 

V;8DSin(6') V 


Id 


Emax is the largest energy electron and Emin = iTie is 
the lowest energy electron that can be emitted in a semi- 
leptonic decay in the rest frame of the heavy meson, and 
Umax is the largest rapidity at which electron measure¬ 
ments are reported (there is an assumption that measure¬ 
ments are restricted symmetrically about y = 0). Note 
that the 9 in the second line of Eq. (Bl) is the usual 
Heaviside step function, and 




(B5) 

(B6) 


The explicit probability functions that were used in 
Eq. (Bl) are (with some numbers truncated slightly for 
visibility, although it should not affect results) 


= {E- eEJ{e^^^ - E){23M23E + ISS.OSdTE^ - 685.7559E^ + 1206.6398E^ 

- 1013.2945^;’^ -h 335.7252£;®)/l.0015 (B7) 

pB^<^{E) = {E- - £;)(0.0015939 - 0.002299E -h 0.069117E2 - 0.089862E3 

-H 0.054773^4 - 0.012189£;®)/0.050435 (B8) 

= -0.003556 -f E(-0.359519 -f E;(27.194146 -h £;(-119.476261 -f E;(244.65394 
-h £;(-292.041775 -f E(218.040963 -h £;(-103.278527 -H E(30.170577 
-h £;(-4.960071 -H £;(0.351141))))))))))(2.65 - £;)/1.004464 (B9) 

pB^J/'l’i^E) = (E(y££^ - 2)(126.32439(£;2 - rn^j/^f + 21.848881(^2 - 
1.728567y££n^ -h 1.721525(E;2 - m^/^) - 64.938832(^2 - 

- 90.585013(^2 - - m2/^E^( 1.728567 - 1.723872E„-H 

32.899797E^ - 63.517521E)^ -h 51.249069E^ - 19.183254E^ -h 2.731110E^)), (BIO) 


where Em = {m% — vr?ji^l2rfiB- See Table I for the 
maximum and minimum allowed energies for each of the 
above P functions. 


TABLE I. Maximum and minimum allowed energies for the 
semi-leptonic electron decay fragments of D and B mesons 
and the J/'i/) mesons from B mesons in the rest frame of 
the decaying heavy meson, rrie = 0.0005 GeV and mj/^ = 
3.096916 GeV. 
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B —>■ e 

B ^ D ^ e 
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me 
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